\documentclass[aps,pre,twocolumn,10pt]{revtex4-1}
%\documentclass[aps,pre,preprint]{revtex4-1}
\usepackage{srcltx}
\usepackage{amsmath}
\newcommand{\Dphi}{{\Delta \phi}}
\newcommand{\Dt}{{\Delta t}}


\begin{document}
 \title{Implementation of the two states model of physiological aging}
 \maketitle

\section{Two states model}

\subsection{Space and time discretization}

The dynamical equation to be numerically implemented is
\begin{equation} \label{difEq}
 \frac{dp_i(\phi,t)}{dt} = - V_i\frac{dp_i(\phi,t)}{d\phi}+\sum_k T_{ik} p_k(\phi,t) .
\end{equation}
with
\begin{equation}
 \sum_i T_{ij} = 0 .
\end{equation}
The number individuals within $\phi_0 \le \phi \le \phi_1$ at state $i$ is
given by
\begin{equation}
 p_i(t) = \int_{\phi_0}^{\phi_1} p_i(\phi,t) d\phi ,
\end{equation}
and the total number of individuals is
\begin{equation}
 p(t) = \sum_i p_i(t) .
\end{equation}
The dynamics of $p(t)$ can be calculated from eq. (\ref{difEq})
\begin{equation}
 \frac{dp(t)}{dt} = \sum_i v_i [p_i(\phi_0,t)-p_i(\phi_1,t)] .
\end{equation}
The above equation defines also the boundary conditions, where $\alpha(t)$
is the density of particles entering in the system at $\phi_0$ and $\alpha'(t)$
is the density of particles leaving the system at $\phi_1$
\begin{align}
\alpha(t) &= \sum_i v_i p_i(\phi_0,t) , \label{bCond0}\\
\alpha'(t) &= \sum_i v_i p_i(\phi_1,t) . \label{bCond1}
\end{align}

We will use a two states model where
\begin{eqnarray}
  V &=& \left[\begin{matrix} 0 \\ \nu \end{matrix}\right] \\
  T &=& \left[\begin{matrix} -k & k \\ k & -k \end{matrix}\right]
\end{eqnarray}

The discretization is done with $N$ points spaced by $\Dphi$
\begin{eqnarray}
t &\rightarrow& n \Dt \\
\phi &\rightarrow& \phi_0 + j \Dphi
\text{ with } \Dphi = \frac{\phi_1-\phi_0}{N-1} \label{dPhi0}\\
\text{or }
\phi &\rightarrow& \phi_0 + j \Dphi + \frac{\Dphi}{2}
\text{ with }  \Dphi = \frac{\phi_1-\phi_0}{N} \label{dPhi1} \\
p_i(\phi,t) &\rightarrow& p_{ij}^n  .
\end{eqnarray}
The discretization of equation (\ref{dPhi0}) has points exactly at the extremes
of the interval [$\phi_0, \phi_1$]. The space is divided in cells of width
$\Dphi$ and the points are at both sides of each cell. Equation (\ref{dPhi1})
is a more conventional scheme, were the points are at the center of each
discretization cell.

\subsection{Upwind differentiation}

In this section the $\phi$ derivative is done upwind, i.e., the derivative at
$p_j$ is computed as
\begin{equation}
 \frac{dp_i(\phi,t)}{d\phi} \rightarrow \frac{p_{ij}^n-p_{i,j-1}^n}{\Dphi} .
\end{equation}
The upwind derivative sometimes give better results when doing forward Euler
differencing. However, we since we performed a Crank-Nicholson like integration
it is not clear if this is the best method to use.

\subsubsection{PDE discretization}

Using the upwind derivative, the discretized form of equation (\ref{difEq}) is
\begin{multline}
 \frac{p_{ij}^{n+1}-p_{ij}^n}{\Dt} =
- V_i \frac{p_{ij}^{n+1}-p_{i,j-1}^{n+1}+p_{ij}^n-p_{i,j-1}^n}{2\Dphi}\\
- k \frac{p_{ij}^{n+1}+p_{ij}^n}{2}+k \frac{p_{1-i,j}^{n+1}+p_{1-i,j}^n}{2} .
\end{multline}
Once the we know the values of $p_{i,j-1}^n$, $p_{ij}^n$ and $p_{i,j-1}^{n+1}$ are know
for both values of $i$, we can calculate the LHS of
\begin{multline}
\left(\frac{1}{\Dt}+\frac{V_i}{2\Dphi}+\frac{k}{2}\right)
p_{ij}^{n+1}-\frac{k}{2}p_{1-i,j}^{n+1} = \\
\frac{V_i}{\Dphi} \frac{p_{i,j-1}^{n+1}+p_{i,j-1}^n}{2} \\
+ \left(\frac{1}{\Dt} - \frac{V_i}{2\Dphi}-\frac{k}{2}\right) p_{ij}^n
+\frac{k}{2}p_{1-i,j}^n .  \label{UWLinSys}
\end{multline}
That system of equations may be written as

\begin{equation}
\begin{split} \label{linSysUpW}
 a_{00} p_{0j}^{n+1} + a_{01} p_{1j}^{n+1} &= b_0\\
 a_{10} p_{0j}^{n+1} + a_{11} p_{1j}^{n+1} &= b_1
\end{split}
\end{equation}
with
\begin{align}
 a_{00} = \frac{2}{k\Dt} + 1 \hspace{1 cm}
 &a_{01} = -1 \\
 a_{10} = -1 \hspace{1 cm}
 &a_{11} = \frac{2}{k\Dt} + \frac{\nu}{k\Dphi} + 1\\
\end{align}
\begin{align}
 b_0 =& \left(\frac{2}{k\Dt} - 1\right) p_{0j}^n + p_{1j}^n\\
 b_1 =& \frac{\nu}{k\Dphi} \left(p_{1,j-1}^{n+1}+p_{1,j-1}^n \right)\nonumber\\
      &+ \left(\frac{2}{k\Dt} - \frac{\nu}{k\Dphi} - 1\right) p_{1j}^n + p_{0j}^n .
      \label{b1}
\end{align}

\subsubsection{Condition for $p_{1j}^{n+1}>0$}

With no inflow in $p_{1j}^n$, i.e., if
\begin{equation}
p_{1,j-1}^n = p_{1,j-1}^{n+1} = p_{0j}^n = p_{0j}^{n+1} = 0 .
\end{equation}
we can rewrite eq. (\ref{UWLinSys}) as
\begin{equation}
\left(\frac{1}{\Dt}+\frac{\nu}{2\Dphi}+\frac{k}{2}\right) p_{1j}^{n+1} =
\left(\frac{1}{\Dt} - \frac{\nu}{2\Dphi}-\frac{k}{2}\right) p_{1j}^n.
\end{equation}
The condition $p_{1j}^{n+1} \ge 0$ must be satisfied even in that extreme case,
implying in
\begin{equation}
\Dt \le \frac{2}{\nu/\Dphi + k} .
\end{equation}

Equivalent analysis of $p_{0j}^{n+1}$ leads to the less restrictive condition
\begin{equation}
\Dt \le \frac{2}{k} .
\end{equation}

\subsubsection{Boundary condition}

The condition (\ref{bCond0}) can be applied to the discretization scheme (\ref{dPhi0}) by doing
\begin{align}
  p_{10}^{n+1} &= \frac{2}{\nu} \alpha(t+\Dt/2) - p_{10}^n\\
  p_{00}^{n+1} &= \frac{1}{a_{00}} \left[\left(\frac{1}{\Dt} - \frac{k}{2}\right) p_{00}^n
    +\frac{k}{\nu}\alpha(t+\Dt/2) \right]
\end{align}
which is equivalent of doing
\begin{align}
 a_{10} &= 0  \\
 a_{11} &= 1 \\
 b_1 &= \frac{2}{\nu} \alpha(t+\Dt/2) - p_{10}^n
\end{align}
However, applying the boundary condition in that way leads to spurious density fluctuation
close to $\phi_0$.

When discretization scheme eq. (\ref{dPhi1}) is used, condition (\ref{bCond1})
can be incorporated by making the substitution
\begin{equation}
p_{1,j-1}^{n+1} + p_{1,j-1}^n \rightarrow \frac{2}{\nu} \alpha(t+\Dt/2) ,
\end{equation}
in equation (\ref{b1}) with $j=0$. The result is that equation
(\ref{linSysUpW}), for $j=$ should use
\begin{align}
 b_1 =& \frac{\nu}{k\Dphi} \frac{2\alpha(t+\Dt/2)}{\nu}\\
      &+ \left(\frac{2}{k\Dt} - \frac{\nu}{k\Dphi}
      - 1\right) p_{1j}^n + p_{0j}^n .
\end{align}

The outflow at the end of the stage, eq (\ref{bCond1}) can be calculated as
\begin{equation}
 \alpha'(t+\Dt/2) = \frac{\nu}{2} \left(p_{1,N-1}^{n+1}+p_{1,N-1}^n \right) ,
\end{equation}
and it can be used as the boundary condition at the beginning of the next stage.

\subsection{Mixed upwind downwind differentiation}

The symmetry properties of eq. (\ref{difEq}) make it more reasonable to use the
upwind differentiation at $t = i\Dt$ and the downwind differentiation at $t =
(i+1) \Dt$. This small change leads to completely different discrete equations.

\subsubsection{PDE discretization}

Using the upwind derivative, the discretized form of equation (\ref{difEq}) is
\begin{multline}
 \frac{p_{ij}^{n+1}-p_{ij}^n}{\Dt} =
- V_i \frac{p_{i,j+1}^{n+1}-p_{ij}^{n+1}+p_{ij}^n-p_{i,j-1}^n}{2\Dphi}\\
- k \frac{p_{ij}^{n+1}+p_{ij}^n}{2}+k \frac{p_{1-i,j}^{n+1}+p_{1-i,j}^n}{2} .
\end{multline}
Once the values of $p_{i,j-1}^n$, $p_{ij}^n$ and $p_{i,j-1}^{n+1}$ are know
for both values of $i$, we can calculate the LHS of
\begin{multline}
\frac{V_i}{2\Dphi} p_{i,j+1}^{n+1} = -
\left(\frac{1}{\Dt}-\frac{V_i}{2\Dphi}+\frac{k}{2}\right) p_{ij}^{n+1}
+\frac{k}{2}p_{1-i,j}^{n+1}  \\
+ \frac{V_i}{2\Dphi} p_{i,j-1}^n + \left(\frac{1}{\Dt} -
\frac{V_i}{2\Dphi}-\frac{k}{2}\right) p_{ij}^n +\frac{k}{2}p_{1-i,j}^n .
\label{sysLinUD}
\end{multline}
The resulting equation for $i=0$ is
\begin{multline}
\left(\frac{1}{\Dt}+\frac{k}{2}\right) p_{0j}^{n+1} = \frac{k}{2}p_{1j}^{n+1} +
\left(\frac{1}{\Dt} - \frac{k}{2}\right) p_{0j}^n +\frac{k}{2}p_{1j}^n ,
\end{multline}
therefore, if we know $p_{1j}^{n+1}$ we can calculate $p_{0j}^{n+1}$. Once
$p_{0j}^{n+1}$ and $p_{1j}^{n+1}$ are known can find the value of $p_{1,j+1}^{n+1}$ by making
$i=1$ in eq. (\ref{sysLinUD})
\begin{multline}
\frac{\nu}{2\Dphi} p_{1,j+1}^{n+1} = -
\left(\frac{1}{\Dt}-\frac{\nu}{2\Dphi}+\frac{k}{2}\right) p_{1j}^{n+1}
+\frac{k}{2}p_{0j}^{n+1}  \\
+ \frac{\nu}{2\Dphi} p_{1,j-1}^n + \left(\frac{1}{\Dt} -
\frac{\nu}{2\Dphi}-\frac{k}{2}\right) p_{1j}^n +\frac{k}{2}p_{1j}^n .
\end{multline}

\subsection{Second-order space accuracy (Crank-Nicolson)}

The second order accurate discretization of the $\phi$ derivative of $P(\phi,t)$ is
\begin{equation}
 \frac{dp_i(\phi,t)}{d\phi} \rightarrow \frac{p_{i,j+1}^n-p_{i,j-1}^n}{2\Dphi} .
\end{equation}
We expect that such a symmetric form will perform better when using a
time symmetric scheme like the Crank-Nicholson. Using that expression
in eq. (\ref{difEq}) with discretization scheme (\ref{dPhi1}) is
equivalent to taking the appropriate average among each pair of neighbors points
to calculate the incoming and outgoing flow in each cell.

\subsubsection{PDE discretization}

Using the second order $\phi$ derivative, the discretized form of equation (\ref{difEq}) is
\begin{multline}
 \frac{p_{ij}^{n+1}-p_{ij}^n}{\Dt} =
- V_i \frac{p_{i,j+1}^{n+1}-p_{i,j-1}^{n+1}+p_{i,j+1}^n-p_{i,j-1}^n}{4\Dphi}\\
- k \frac{p_{ij}^{n+1}+p_{ij}^n}{2}+k \frac{p_{1-i,j}^{n+1}+p_{1-i,j}^n}{2} .
\end{multline}
Once the we know the values of $p_{i,j-1}^n$, $p_{ij}^n$ and $p_{i,j-1}^{n+1}$ are know
for both values of $i$, we can calculate the LHS of
\begin{multline}
-\frac{V_i}{2k\Dphi} p_{i,j-1}^{n+1} +  \left(\frac{2}{k\Dt}+1\right) p_{ij}^{n+1} +\frac{V_i}{2k\Dphi} p_{i,j+1}^{n+1} \\
- p_{1-i,j}^{n+1} =
\frac{V_i}{2k\Dphi} (p_{i,j-1}^n-p_{i,j+1}^n) \\
+ \left(\frac{2}{k\Dt} -1\right) p_{ij}^n +p_{1-i,j}^n. \label{CNEquation}
\end{multline}
That system of equations may be written as

\begin{equation}
 a  P_j^{n+1} = b P_j^n
\end{equation}
with
\begin{equation}
P_j^n = \left[\begin{array}{l}
  p_{0,j-1}^n \\ p_{1,j-1}^n \\ p_{0j}^n \\ p_{1j}^n \\ p_{0,j+1}^n \\ p_{1,j+1}^n
\end{array}\right] ,
\end{equation}

\begin{equation}
a = \left[\begin{array}{cccccc}
     0 & 0 & \frac{2}{k\Dt}+1 & -1 & 0 & 0  \\
     0 & -\frac{\nu}{2k\Dphi} & -1 & \frac{2}{k\Dt}+1 & 0 & \frac{\nu}{2k\Dphi}
    \end{array}\right] \label{CNa}
\end{equation}

\begin{equation}
b = \left[\begin{array}{cccccc}
     0 & 0 & \frac{2}{k\Dt}-1 & 1 & 0 & 0  \\
     0 & \frac{\nu}{2k\Dphi} & 1 & \frac{2}{k\Dt}-1 & 0 & -\frac{\nu}{2k\Dphi}
    \end{array}\right] . \label{CNb}
\end{equation}

The whole problem can be written as the following linear system involving
$2N\times 2N$ matrices and $2N$-vectors:
\begin{equation}
 A  P^{n+1} = B P^n \label{CNLinSys}
\end{equation}
with
\begin{equation}
P^n = \left[\begin{array}{l}
  p_{0,0}^n \\ p_{1,0}^n \\ p_{01}^n \\ p_{11}^n \\ \vdots \\ p_{0,N-1}^n \\ p_{1,N-1}^n
\end{array}\right] ,
\end{equation}
\begin{equation}
A = \left[\begin{array}{cccccccccc}
     a_{02} & a_{03} & a_{04} & a_{05} & 0 & 0 & 0 & 0 &  & \\
     a_{12} & a_{13} & a_{14} & a_{15} & 0 & 0 & 0 & 0 &  & \\
     a_{00} & a_{01} & a_{02} & a_{03} & a_{04} & a_{05} & 0 & 0 &  & \\
     a_{10} & a_{11} & a_{12} & a_{13} & a_{14} & a_{15} & 0 & 0 &  & \\
     0 & 0 & a_{00} & a_{01} & a_{02} & a_{03} & a_{04} & a_{05} &  & \\
     0 & 0 & a_{10} & a_{11} & a_{12} & a_{13} & a_{14} & a_{15} &  & \\
       &   &        & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \\
     & & & & a_{00} & a_{01} & a_{02} & a_{03} & a_{04} & a_{05}\\
     & & & & a_{10} & a_{11} & a_{12} & a_{13} & a_{14} & a_{15}\\
     & & & & 0 & 0 & a_{00} & a_{01} & a_{02} & a_{03}\\
     & & & & 0 & 0 & a_{10} & a_{11} & a_{12} & a_{13}\\
     \end{array}\right] ,
\end{equation}
and similarly for $B$.

Since $a_{00}=a_{10}=0$, matrix $A$ have, at most, two non-zero elements to the
left of the diagonal and three non-zero elements to the right of the diagonal.
The computational time necessary for solving such band diagonal matrices goes
as $N$, when specialized routines are used.

\subsubsection{Boundary condition}

The boundary condition (\ref{bCond0}) can be implemented by replacing
\begin{equation}
p_{1,j-1}^{n+1} + p_{1,j-1}^n \rightarrow \frac{2}{\nu} \alpha(t+\Dt/2) ,
\end{equation}
in eq. (\ref{CNEquation}) when $j=0$. The only difference on matrices $a$ e $b$
of equations (\ref{CNa}) and (\ref{CNb}) due to that repalcement
 is that the second element of the
second line will be zero when $j=0$. For that value of $j$, instead of equation
(\ref{CNLinSys}), we must use
\begin{equation}
aP_0^{n+1} = b P^n_0 + \left[\begin{array}{c} 0 \\ 1\end{array}\right]
\frac{1}{k\Dphi} \alpha(t+\Dt/2) .
\end{equation}

At the right boundary we must use the upwind $\phi$ derivative, eq.
(\ref{UWLinSys}), instead of eq. (\ref{CNEquation}). Therefore, when $j=N-1$,
the matrices to be used are
\begin{equation}
a = \left[\begin{array}{cccccc}
     0 & 0 & \frac{2}{k\Dt}+1 & -1 & 0 & 0  \\
     0 & -\frac{\nu}{k\Dphi} & -1 & \frac{2}{k\Dt}+\frac{\nu}{k\Dphi}+1 & 0 & 0
    \end{array}\right]
\end{equation}
\begin{equation}
b = \left[\begin{array}{cccccc}
     0 & 0 & \frac{2}{k\Dt}-1 & 1 & 0 & 0  \\
     0 & \frac{\nu}{k\Dphi} & 1 & \frac{2}{k\Dt}-\frac{\nu}{k\Dphi}-1 & 0 & 0
    \end{array}\right] .
\end{equation}
The outflow of particles seems to be
\begin{equation}
\alpha'(t+\Dt) = \frac{\nu}{2}\left(p_{1,N-1}^{n+1}+p_{1,N-1}^n\right) . \label{wrongAlpha}
\end{equation}

However, particle conservation imposes that the terms coming from $dp_1(\phi,t)/d\phi$
must cancel among different lines, except by the rightest terms, which contribute
to the outgoing flow. That contribution can be measured by
\begin{equation}
\alpha'(t+\Dt/2) = -\frac{\Dphi}{\Dt} \sum_j \Delta p_{1,j}^n , \label{alphaDelta}
\end{equation}
where $\Delta p_{1,j}^n$ is the variation due to discretization of  $dp_1(\phi,t)/d\phi$.
The unbalanced term originate in line $N-2$ is
\begin{equation}
\Delta p_{1,N-2}^n = -\frac{\nu\Dt}{4\Dphi} \left(p_{1,N-1}^{n+1}+p_{1,N-1}^n\right),
\end{equation}
and the unbalanced term originated in line $N-1$ are
\begin{equation}
\begin{split}
\Delta p_{1,N-1}^n =& -\frac{\nu\Dt}{2\Dphi}\left(p_{1,N-1}^{n+1}+p_{1,N-1}^n\right)\\
 &+\frac{\nu\Dt}{2\Dphi}\left(p_{1,N-2}^{n+1}+p_{1,N-2}^n\right)
\end{split}
\end{equation}
Putting these contributions together in eq. (\ref{alphaDelta}) we arrive to
\begin{equation}
\begin{split}
\alpha'(t+\Dt) =& \frac{3\nu}{4}\left(p_{1,N-1}^{n+1}+p_{1,N-1}^n\right)  \\
&- \frac{\nu}{4}\left(p_{1,N-2}^{n+1}+p_{1,N-2}^n\right),
\end{split}
\end{equation}
which should be used instead of eq. (\ref{wrongAlpha}).

\section{Maximum Likelihood}

\subsection{Individual Completion Time}

The times $t_0$ and $t_1$, at which a given individual started and finished
the phase, are known to be in intervals $[t_{00}, t_{01}]$ and $[t_{00}, t_{01}]$.
The distribution of these variables are assumed to be
\begin{equation}
p(t_i) =
 \begin{cases}
    0  & \text{if } t<t_{i0} \\
    (t_{i1}-t_{i0})^ {-1}   & \text{else if } t<t_{i1}\\
    0  & \text{otherwise}\\
 \end{cases} .
\end{equation}
The distribution of the completion time
\begin{equation}\
	\tau = t_1 - t_0
\end{equation}

\begin{align}
	p(\tau) & = \int p_0(t) p_1(t+\tau) dt \\
	 = &  \begin{cases}
    	0  & \text{if } \tau<t_{10}-t_{01} \text{ else} \\
    	\frac{\min(t_{11}-\tau,t_{01})-\max(t_{10}-\tau,t_{00})}{(t_{01}-t_{00})(t_{11}-t_{10})}
     	& \text{if } \tau<t_{11}-t_{00}\\
    	0  & \text{otherwise }
	 \end{cases}
\end{align}
which may be written as
\begin{widetext}
\begin{align}
p(\tau) = &  \frac{1}{(t_{01}-t_{00})(t_{11}-t_{10})} \times \begin{cases}
    	0  & \text{if } \tau<t_{10}-t_{01} \text{ else} \\
    	\tau-t_{10}+t_{01}& \text{if } \tau<\min(t_{10}-t_{00},t_{11}-t_{01}) \text{ else} \\
    	\min(t_{01}-t_{00},t_{11}-t_{10}) & \text{if } \tau<\max(t_{10}-t_{00},t_{11}-t_{01}) \text{ else} \\
    	t_{11}-t_{00}-\tau & \text{if } \tau<t_{11}-t_{00}\\
    0  & \text{otherwise }  .
 \end{cases}
\end{align}
\end{widetext}

The above may be rewritten is we definde
\begin{equation}
\begin{split}
\tau_1 &= t_{10}-t_{01}\\
\tau_2 &= \min(t_{01}-t_{00}, t_{11}-t_{10})\\
\tau_3 &= t_{11}-t_{00}
\end{split}
\end{equation}
%\begin{widetext}
\begin{align}
p(\tau) = &  \frac{1}{\tau_3-\tau_1-\tau_2} \times \begin{cases}
    	0  & \text{if } \tau<\tau_1 \text{ else} \\
    	\frac{\tau-\tau_1}{\tau_2} & \text{if } \tau<\tau_1+\tau_2 \text{ else} \\
    	1 & \text{if } \tau<\tau_3-\tau_2 \text{ else} \\
    	\frac{\tau_3-\tau}{\tau_2} & \text{if } \tau<\tau_3\\
    0  & \text{otherwise }  .
 \end{cases}
\end{align}
%\end{widetext}



\end{document}
